#include "stdafx.h"
#include "basic_data.h"

using namespace user_constant;

double get_momentum_energy(double energy_MeV){return sqrt(energy_MeV*energy_MeV - me_MeV*me_MeV);}
double get_momentum_beta(double beta){return beta*sqrt(1/(1-(beta*beta)))*(me_MeV/c);}
double get_momentum_gamma(double gamma){return get_momentum_energy(gamma*me_MeV);}
double get_momentum_kinetic(double w){return get_momentum_energy(w+me_MeV);}

double get_kinetic_momentum(double cp_MeV){return sqrt((cp_MeV*cp_MeV - me_MeV*me_MeV) - me_MeV);}
double get_kinetic_betagamma(double betagamma){return get_kinetic_momentum(betagamma*me_MeV/c);}

void convert_to_phase(phase_space_point& output,trace_space_point& input,double z0)
{
	output.x = input.x;
	output.y = input.y;
	output.z = input.phi/360*(c/f0) - z0;
	
//	const double p = get_momentum_kinetic(input.w)/(me_MeV/c);
	const double gamma = (input.w + me_MeV)/me_MeV;
	const double beta = sqrt(1-1/(gamma*gamma));
	const double p = beta*gamma;

	output.px = p*input.xp;
	output.py = p*input.yp;
	output.pz = sqrt(p*p - output.px*output.px - output.py*output.py);

}

void convert_to_phase(phase_space_dist& output,trace_space_dist& input,double z0)
{
	output.clear();
	output.set(input.size);

	for(int i=0;i<input.size;i++)
	{
		convert_to_phase(output.dataset[i],input.dataset[i],z0);
	}
}

void convert_to_trace(trace_space_dist& output,phase_space_dist& input)
{
	output.clear();
	output.set(input.size);

	double momentum = 0;

	for(int i=0;i<input.size;i++)
	{
		momentum = input.dataset[i].p();

		output.dataset[i].x = input.dataset[i].x;
		output.dataset[i].y = input.dataset[i].y;
		output.dataset[i].phi = (input.dataset[i].z - input.dataset[0].z)/(c/f0)*360;

		output.dataset[i].xp = input.dataset[i].px/momentum;
		output.dataset[i].yp = input.dataset[i].py/momentum;
		output.dataset[i].w = get_kinetic_betagamma(momentum);

	}
}